Inhomogeneous vortex matter 
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We present a generalization of the continuum theory of vortex matter for non-uniform superfluid 
density. This theory explains the striking regularity of vortex lattices observed in Bose-Einstein 
condensates, and predicts the frequencies of long-wavelength lattice excitations. 



Dense lattices of quantized vortices in rotating Bose- 
Einstein condensates (BECs) jj], ^ are strikingly more 
regular than finite vortex arrays in homogeneous super- 
fluid § (see Fig. 1), even though BEC densities vary 
greatly over the sample. This Letter generalizes the 
Feynman-Tkachenko Q ^) continuum theory of 'vortex 
matter' to cases in which the condensate density varies 
slowly on the scale of the lattice spacing. This theory 
explains the lattices' surprising regularity, and find pro- 
nounced effects of nonuniform density on lattice excita- 
tions. 
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FIG. 1: (a) Static lattice according to Eq. (p^|) , translated and 
rotated to match (b) experimental data courtesy of J.R. Abo- 
Shaeer. Compare (c) vortex array in constant p (2172 from 
Fig. 5 of Ref. §). 



We consider a two-dimensional regular array of vor- 
tices, whether realized in a very oblate BEC, or as parallel 
vortex lines in a prolate one. Denoting the lattice length 
scale by b and introducing dimensionless complex co- 
ordinates bz = x + iy, a regular lattice has the positions 
zjk of parallel vortex lines given by Zjk — Zj k = k + rj 
for r = n + iT2 (t2 > and r, real). Much is known 
about vortex lattices in superfluids of constant density 
p, simply from incompressible hydrodynamics; and it is 
all simplified by using dimensionless variables, express- 
ing time, velocity, and energy in lattice units Mb 2 /h, 
H/(Mb), h 2 /(Mb 2 ) respectively, for M the mass of the 
particles composing the superfluid. The regular trian- 
gular case t = (1 + iV3)/2 is the ground state 0] of a 
sample rotating at dimensionless rate Q = tt/t2. The ir- 
rotational velocity field v = v x + iv y consists of a fine field 
Vf, which is periodic on the lattice scale, plus a coarse 



field v c obeying Feynman's criterion 

d z v c - d s v c = 2-Kipv (1) 

where for a regular lattice the vortex density pv is 1/t2. 
(In our complex notation 2d z A = 6(V • A + iV X A) for 
any A. We assume counter-clockwise rotation 
long- wavelength excitations of the lattice, zjk(t) 

D (^ z jk>Zj k ; tj obey the wave equations Q 

id t {d z D~d- z D) = 2fl(d z D + d z D) 
1 



And 



id t (8 Z D + d- z D) 



-d- zz (d z D - d- z D) 



(2) 
(3) 
be- 



(Extensions to three-dimensional vortex matter | 
come considerably more complicated.) 

Although current dilute gaseous BECs are compress- 
ible fluids governed by the Gross-Pitaevskii equation 
(GPE) ||, much of BEC vortex physics can be cast into 
a simpler hydrodynamic form. In our dimensionless vari- 
ables, the GPE in a frame rotating about the co-ordinate 
origin may be written 



d t p 
d t v 



fld^p - [d z (pv) + d z (pv)} 
d z \v-iflz\ 2 + 2V + 2gp- 



dzz^fp 



VP 



(4) 
(5) 



V = V t 



n 2 2 

trap 2~ l^l 



where v is still the lab-frame velocity, <f> is the polar an- 
gular co-ordinate, Vt rap (z, z) is the trap potential, and g 
is the dimensionless 2D coupling, determined by atomic 
and trap parameters. Except very near vortex cores, 
the rotating-frame velocity v — iflz is of order unity; 
and p~ x l 2 d zz ^fp is no larger, except in cores and sam- 
ple edges, which may be treated separately as bound- 
ary layers. In current experiments the healing length 
£ = (gp) 1 , which sets the vortex core size, is every- 
where else much smaller than b. 

So, outside vortex cores, the leading order results in 
the co-rotating frame are 

V 

p = const. (6) 



.9 



fld^p = d z (pv) + d z (pv) . 



(7) 
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Thus p depends on v only through the centrifugal mod- 
ification of the trap potential. The velocity field, v is 
determined by (]?]), plus the condition of irrotationality 
except at the quantized vortex cores: 



d z v - d z v 



1m ( z 

jk 



Zjk) 



(8) 



This, with V • v = instead of (Q), is the starting point 
for Ref. [||. So for constant V, and hence uniform p, the 
results of the incompressible case also apply to BECs. 

If p varies slowly in space, will not inhomogeneous ef- 
fects be small? Not obviously: like p itself, the lattice 
shape might vary slowly, but change greatly over the 
whole sample. Indeed, even for constant p in the ground 
state of a finite vortex array, Campbell and Ziff || found 
gradual but cumulatively large distortions. But there is a 
basic problem in extending their analysis to non-uniform 
P- 

Investigations of homogeneous vortex matter have gen- 
erally relied on the exact single-vortex solution v — 
vi(z — Zjk) to (0) and (||) for constant p, which one 
may simply sum over the vortex labels j, fc, because (Q) 
and @ are linear in v. For non-constant p the famil- 
iar vi(z) = i/z satisfies @ but not (fj]), and so we do 
not have the exact single-vortex solution. Pcrturbative 
approximations about v\{z) = i/z as an ansatz break 
down at distances from the core beyond the length scale 
of the density variation JsJ . So the few- vortex problem in 
inhomogeneous BECs becomes analytically intractable. 
For sufficiently dense lattices, however, inhomogeneous 
vortex matter yields to a different approach. 

By a 'dense' vortex lattice, we mean p — p(ez,ez) for 
small e. (For a round harmonic trap with Thomas-Fermi 
radius R, e = b/R, giving s ~ 0.1 in current experi- 
ments.) We can therefore perturb in e; but to distinguish 
smallness from slowness, we must use multiple scale anal- 
ysis (MSA) |IC[ |. This formalism produces coarse-scale 
equations of motion for D, from which all lattice-scale 
physics has been eliminated, in the same sense that high- 
frequencies are eliminated by adiabatic methods. These 
will be our generalizations of (Q) and (||). 

The application of MSA leads to a rather involved 
derivation the details of which will be reported elsewhere; 
here we outline its steps and report its conclusions. We 
begin by satisfying (0) identically by defining 

v = iQz + -d z (pF) (9) 
P 

for real F. We use vortex-centered co-ordinates: 

z = z' + D(z',z',t), (10) 
so that z\ h = z\ regardless of D. We do this so that 



in the z' co-ordinates we always have a regular lattice, 
whose symmetries we can exploit, even though the phys- 
ical lattice may be distorted by excitations, or by a 



static D field induced by inhomogeneous p. In the non- 
Cartesian z' co-ordinates, Eqn. M) becomes 



Re 



dz pF I = 7T 2^ o {z -zf, 



jk 



Q l + d*,D + d s >D- 



d zJ D d z >D 
d z ,D d~>D 



This shows explicitly that only gradients in D affect F. 

MSA then embeds the physical z'-plane in a fictitious 
4-space of complex co-ordinates f , Z, as the subspace 
(C,Z) — (z',ez'). This provides d z > — » 8q + edz, etc., 
and proceeding perturbatively in e, we are able to write 
explicit solutions (in terms of rapidly converging series) 
for the ^-dependence of F at every order (we need to 
go to third). As usual with MSA, the 'gauge freedom' 
in how functions depend explicitly on the two extra di- 
mensions is used to remove solutions growing secularly 
with C, by constraining the purely Z-dependent part of 
F (which we denote by F c ). Once we restrict back to 
physical two-space by setting Z — ez', C = z' ', and re- 
turn from vortex- fixed z' to Cartesian z, we recognize 
the constraint on F c as just what is needed to maintain 
Feynman's condition (|l]), for py as perturbed by D: 



0: 



dz (pF c ) 



P 



d _WV = - 2 n (o z d + d- z D) , (ii) 



where we drop determinant terms quadratic in D (be- 
cause we will have D order e or smaller). 

Having solved for F, and hence v, in terms of explicit 
lattice-periodic functions and v c , we know the local fluid 
velocity near each vortex. This fixes the instantaneous 
vortex translational velocity field D; but the fixing is not 
trivial. Since the hydrodynamic approximation to the 
GPE breaks down within \z — Zjk\ ~ £/b, in these small 
regions we must solve the time-dependent GPE using a 
different perturbation theory, based on Taylor-expanding 
V about Zjk- Matching the hydrodynamic and core solu- 
tions smoothly together (see Refs. j^, 11 ) finally yields, 
to leading order in e, 

i P b = x - [d- z ( P d z D + P d- Z D) - d z ( P d z D)] - d- z ( p f c ) 



(d- zP y 



In— + 1.17 &p- 0.20$ 
2tt£ J p 



(12) 



which is expressed in the co-rotating frame. 

Eqns. ([ll]) and ( fl2|) are our main results. The numer- 
ical co-efficients in (|12|) include some numerically evalu- 
ated contributions from the nonlinear core regions (com- 
pare with p!l|), and also functions of r, generally re- 
lated to ^-functions, evaluated for the triangular case 
t = (1 + iv3)/2. For general r (i.e. for lattices other 
than the regular triangular) , (|l^) would have several ad- 
ditional terms, such as B(r)d zz D, where B(t) is another 
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rapidly converging series. Modular covariance (a general 
type of lattice symmetry) of the extra co-efficients like 
B(t) constrains them to vanish when r = 1 +| 

MSA implies that the lattice scale b is to the vortex 
matter equations much as the healing length £ is to the 
hydrodynamic equations that underly them. Thus Eqns. 
( pi] ) and Jl^ ) should be accurate except at distances of 
order b or less from of the edge of a vortex array. But is 
this claim really compatible with the results of Campbell 
and Ziff ||) for finite vortex arrays in infinite homoge- 
neous superfluid? Setting p constant, to leading order 
in e the pair @ dll|) reduces to (§) and @), and set- 
ting D — > Dq = 0, we find a wealth of solutions to these 
fourth order equations: 



Do 



E 

rn — 1 



{a m z m - x + (m + l)b m z m z - b* m z m+1 ) (13) 



for arbitrary complex constants a m , b m . A priori it is un- 
clear what boundary condition D should respect at the 
array edge; however if we assume that the edge should be 
a uniform circle of vortices, fitting leads to a unique com- 
bination of multipolar distortions with m = 6,12,18.... 
Figure 2 shows that stopping at only m = 12 gives quite 
good agreement with the first excited state for 217 vor- 
tices found in Ref. ||. (The ground state differs only in 
the outermost ring.) 



a) * * * I * * * 

/ + + + * + + + 



b) 



Using (|j|) and (|lj), the v c associated with this Dq is 



VcO 



n - 



e 2 fln 



27rg 



1.17] ' 



1-e 2 z 



(15) 



This purely azimuthal flow is slightly less than rigid body 
rotation at Q, and the magnitude of the backflow in- 
creases with radius. This agrees qualitatively with the 
numerical results shown in Fig. 5 of Ref. fTs^| . 

Any static distortions forced by boundary conditions, 
like those of Fig. (2), would appear as zero- frequency 
modes among the collective excitations; so to these we 
now turn. We write D = Do + dg(P + iQ), without loss 
of generality, for real P, Q. Then Re{<9 z [@/p] } and @ 
yield 



d- zz Q = -2Sld- zz P x [1 + 0(e 2 )] 



(16) 



Since only D is physical, any terms in P annihilated by 
dzz are of form f(z) + f*(z) and so can be absorbed in 
Q as i[f*(z) — f[z)]. And since one can easily show that 
Laplacian-free terms in Q must be time-independent to 
leading order in e 2 , we can set P = —^Q- 

Introducing polar co-ordinates re 1 ^ = ez, so that p = 
po(l—r 2 ), we can write Q = g mri (r) cos (mcfi — u> mn t — a) 
(for arbitrary constant a and angular and radial quantum 
numbers m and n). Considering first the case m = 0, 
where d^Fc = 0, the imaginary part of e~^x (|l^) implies 



♦ + * * * ♦ 
+ + + + + + 
+ + + + + 



+ + + + + + + 

+ + + 4- + + 
+ +♦ + + + + 

+ + + + * + «- 
* «* + * * * * i 



* • 



FIG. 2: Distortions of a finite vortex array in infinite homoge- 
neous superfluid. Fig. 2a) shows the combination of m — 0, 6, 
and 12 solutions that makes the outer ring most circular and 
evenly spaced, while 2b) is a) overlaid on array 2172 of Fig. 
5 of Ref. §. 

If we set p = po(l — s 2 \z\ 2 ) for a BEC in a round 
harmonic trap, though, the only distortion forced on us 
by the inhomogeneity is the very mild 



Dn = 



V3 



4tt 1 



In— + 1.17) +G(e 4 ) (14) 



This scarcely visible radial shift of each vortex is shown 
in Fig. 1, for e = 1/6 = £(0)/&. For vortices very close to 
the TF surface where formally £ — > oo, ( fl4| ) spuriously 
predicts large inward displacements. (Four such vortices 
have been excised from Fig. la.) Apart from this fail- 
ure in the lattice-edge boundary layer of thickness b, the 
accord with experiment is excellent. 



lHl drqan 



2r" 



(1-r 2 ) 



[(I-. 2 ) (d r - 



d r qon 



(17) 

plus order e. Frobenius analysis shows that only one so- 
lution to this second order equation for d r qo n is finite at 
r — > 0, and imposing finiteness at r — > 1 as well fixes 
the discrete spectrum. (Modified behavior within the 
boundary layer r > 1 — e will be able to reconcile our 
coarse-scale D and v c with microscopic boundary condi- 
tions, as long as our functions remain regular; compare 
the hydrodynamic derivation of collective modes in the 
vortex-free BEC Numerical search yields 



eft{0, 1.43, 2.32, 3.18, 4.03,...}. 



(18) 



The zero mode is global rotation of the lattice, goo = r 2 ■ 

These results are unsurprising for Tkachenko waves in 
a finite cylindrical system; but for m 7^ 0, the dynamics is 
very different. To eliminate F c we must take Im[(9 2 (p"2|)], 
and the non-constant p leaves a first order time derivative 
on the left side, implying u) mn of order e 2 instead of e. 
Our leading order equation is thus first order in t but 
fourth in r: 



16muj ri 



1 r / q 

~Qmn — [(^T? 



r 1 d r - m 2 r 2 ) - Ard r ] 



^2 -2\ 

m r ) q„ 



(19) 
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This equation is quite singular (Frobenius analysis shows 
that of the four solutions only one is finite at both 
r = 0,1); but its differential operator is Hermitian and 
self-dual, and its regular solution has a rapidly con- 
verging Frobenius series in r. From (|TJ) we see that 
the radial component of v c blows up at r = 1 unless 
{d rr — r~ 1 d r — m 2 r~ 2 ) q mn (l) = 0. This boundary con- 
dition on the regular solutions q mn (r) fixes the discrete 
spectrum, which may be found by summing Frobenius 
series numerically: 



LU mn 


m = 1 


2 3 4 5 6 


n = 





0.365 0.900 1.60 2.46 3.49 


1 


2.93 


5.02 7.53 10.5 13.8 17.5 


2 


32.0 


31.3 35.6 41.6 48.8 56.8 


3 


130. 


105. 106. 113. 124. 137. 



and so on. The translational zero mode is q±o = r. These 
eigenvalues of a fourth order equation increase rapidly 
with radial quantum number n, and as the co-efficients 
reach order e" 1 (19) becomes invalid, and WKB-like 
Tkachenko waves will emerge instead. 

The only zero-frequency solution satisfying the bound- 
ary conditions is the rigid translation; and so ( |l4| ) is the 
full static distortion. This reflects the fact that p is so 
small at the edges of the lattice that no boundary ener- 
gies are large enough to influence the bulk lattice. 

Only positive m have been reported, because oj_ mi „ = 
—Lo m n and replacing m — > — m leaves our ansatz for Q un- 
changed. What this means is that the nonuniform p has 
drastically split the degeneracy of the two modes that 
would, for constant p, be proportional to e ±m ^, with 
frequency of order e. The linear combination propor- 
tional to cos(m(j)—uj mn t) propagates much more slowly, in 
the co-rotating frame; evidently the orthogonal combina- 
tion, which we have not examined, propagates much more 
quickly. (The distortion patterns rotate about the origin; 
the vortices follow elliptical orbits about their equilib- 
rium positions: see Fig. 3.) Thus the lowest frequency 
lattice modes have a first order dynamics, and only half 
as many distinct modes as for constant p. 




FIG. 3: Vortex lattice excitations in a round harmonic trap. 
The grey 'trails' indicate vortex motion. Left: m,n — 2,0; 
compare Figure 3 b) of Ref. Q|. Right: m, n = 0, 1, in which 
motion is almost purely angular, and much faster. 



Finally, note that for the quadrupole mode u>2o ^ 
(the zero eigenvalue solution 920 = 7,2 does not satisfy 
the boundary condition). Since it is this mode which 
would distort the equilateral triangular lattice into the 
moderately different regular lattices that are also dynam- 
ically stable on short wavelengths [EJ, we conclude that 
although stable in bulk those lattices are frustrated in 
the finite system, and cannot even be stationary with- 
out slow but cumulatively large distortions. Reviewing 
our calculations in this context, it is clear that the only 
reason the equilateral triangular lattice does not suffer 
a similar fate is the vanishing, due to lattice symmetry, 
of several awkward terms from (UJ) . Once the threat of 
cumulative distortion is lifted, it is not surprising that 
merely local distortion is of order e 2 . So the regularity 
of the observed BEC vortex lattices is ultimately due to 
their triangular structure. An engineer would attribute 
this to triangular rigidity, and a mathematician to the 
fact that the triangular lattice is the Z3 fixed point of 
the modular group. Physicists are entitled to arbitrary 
linear combinations of the two explanations. 
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